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Suinnary 

This  paper  presents  initial  results  on  spatially 
variant  recursive  estimation  of  images.  Parameters 
of  block-wise  constant  recursive  model  are  identif- 
ied on  noise  free  data.  The  models  are  then  used 
to  design  reduced  update  Kalman  filters  which  are 
applied  to  noisy  data.  The  results  are  presented 
and  discussed. 

^ Introduction 

Digital  processing  of  images  has  in  recent  years 
become  both  economic  and  practical.  Most  sophisti- 
cated image  processing  is  performed  off  line  on 
large  machines  because  of  the  large  memory  and  com- 
putational requirements  of  the  presently  used  non- 
recursive methods.  In  particular  for  the  image 
estimation  problem,  classical  nonrecursive  techni- 
ques involve  operations  with  large  matrices  and 
their  Inverses  and  are  hence  not  suitable  for  real- 
time applications  which  might  include: 

1)  Restoration  of  noisy  images  after  reception 
on  a low  power  transmission  link. 

2)  Pictures  arising  from  low  light  level  imaging 
where  back  sensor  noise  significantly  con- 
tributes to  the  output  signal. 

3)  Reception  of  a decoded  DPCM  image  which  re- 
sults from  a maximum- likelihood  decoding 
technique. 

4)  Processing  of  non-image  two-dimensional  (2-D) 
data  for  noise  reduction  prior  to  display  in 
image  format. 

Previous  efforts  towards  the  development  of  re- 
cursive 2-D  filters  have  resulted  in  algorithms 
which  for  the  most  part  take  advantage  of  one-dimen- 
sional approximations  or  require  a state  vector  with 
an  exceptionaly  large  number  of  components  (>100). 
Furthermore  these  algorithms  do  not  take  into  ac- 
count the  non-stationarity  of  the  image  being  pro- 
cessed. Thus  a constant  coefficient  two-dimensional 
model  has  in  general  been  apriori  assumed  for  the 
entire  picture. 

Consequently,  in  view  of  the  need  for  developing 
a two-dimensional  recursive  filter  suitable  for  pro- 
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cessing  non-stationary  images,  we  have  combined  a 
least  squares  parameter  Identification  procedure 
with  the  previously  developed  2-0  reduced  update 
Kalman  filter.' 

More  precisely  we  consider  image  models  as  given 
by  a half-plane  Markovian  model^,  which  is  ob- 
served in  the  presence  of  an  additive  white  Gaus- 
sian noise  field.  Given  a set  of  measurements  and 
values  for  the  model  coefficients,  a Kalman  type 
estimator  can  be  designed  for  estimating  the  image 
pixel  intensities.'  Computation  can  be  greatly 
reduced  by  limiting  the  update  process  at  point 
(m,n)  to  only  those  elements  of  the  image  which  are 
directly  coupled  to  point  (m,n)  via  the  Markovian 
model.  If  however,  the  image  is  non-stationary 
then  the  coefficients  will  vary  as  a function  of 
the  coordinates  (m,n).  In  this  case  it  is  desirable 
that  on-line  identification  of  the  coefficients  be 
performed  in  conjunction  with  the  state  estimation. 

To  this  effect,  least  squares  identification  was 
used  to  find  estimates  of  the  coefficients  over 
blocks  of  the  image.  These  estimates  were  then 
used  by  the  reduced  update  Kalman  filter  for  pro- 
cessing those  image  points  contained  within  the 
pertinent  block.  Results  indicate  that  this  pro- 
cedure significantly  improves  on  the  results  ob- 
tained using  a spatially  invariant  model. 

Reduced  Update  Kalman  Filter 

In  one  dimension,  the  Kalman  filter  offers  an 
attractive  solution  to  the  linear  filtering  and 
prediction  problem.  The  extension  of  one  dimen- 
sional Kalman  filtering  to  two  dimensions  requires 
not  only  a suitable  2-D  recursive  model  but  also 
an  enormous  amount  of  data  storage  and  transfer 
due  to  the  high  dimension  of  the  resulting  state 
vector.  Hence  a straightforward  extension  is  of 
limited  success,  and  thus  it  becomes  desireable  to 
consider  computationally  effective  approximations. 
Here  we  review  one  such  approximation,  the  2-D 
reduced  update  Kalman  filter  as  presented  in  [1]. 

To  Illustrate  this  approach,  consider  the  scan- 
ning of  a discrete  2-D  field  on  an  NXN  regularly 
spaced  lattice.  Since  the  scanning  operation  does 
not  qualitatively  affect  the  results,  we  assume  a 
raster  scan. 

We  now  consider  a signal  model  which  is  Markov- 
ian and  given  by  a non-symmetric  half  plane  (NSHP) 
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recursive  model. 2 

s(m,n)  • I C|^^s(m-k,n-l)  + w(m,n)  (1) 

where  w{ni,^ 
is  an  NSHP, 

this  model  is  {MxM)th  order.  The  observation  model 
is 

r(m,n)  • s(m,n)  + v(m,n)  (2) 

where  v(m,n]  is  a white  Gaussian  source.  Using  the 
scanning  operation  we  transform  the  2-D  problem  in- 
to an  equivalent  1-D  problem.  Define  a state  vec- 
tor of  m(n+l)  components , 

s(m,n)  - [s(m.n),  s(m-l,  n) ,. . . ,s(l  ,n) ; s(n,n-l). 


~ ..$(1 ,n-l );... ;s(n,  n-fflj,...,s{m-m,  n-mj]T 

then  (1)  and  (2)  can  be  put  into  the  form, 

s.(m,n)  »££(m-1,n)  +w(m,n),  (3) 

r(m,n)  • H ^(m.n)  + v(m,n)  (4) 


Thus,  we  could  inmediately  write  down  the  Kalman 
equations'  with  the  above  interpretation  of  the  s 
vector.  The  difficulty  with  these  equations  is  the 
amount  of  computation  and  memory  requirements  asso- 
ciated with  them.  By  limiting  the  update  process 
to  only  those  elements  "near"  the  'present'  point, 
the  computation  can  be  greatly  reduced.  The  re-  • 
suiting  reduces  update  Kalman  filter  equations  can 
be  written  in  scalar  form  as  given  below.  For  de- 
tails see  [1].  In  these  equations,  the  superscript 
Indicates  the  step  in  the  filtering,  while  argu- 
ments represent  the  position  of  the  data  on  the  nxn 
grid. 


) is  a white  Gaussian  noise  field  and^^ 
i.e.  {m>0,  n>0>  \J  {m<0,  n<0}.  Ue  assimle 


State  Prediction  and  update: 

«{m,n)  «(m-l,n) 

Sulm.n)  « 7 c.^  s (m-k,n-i)  (5) 

s^"'’'’>(1,j)  • s^"'’'’>(i,j)  + 

K('"*'’Wi.n-j)-[r(m,n)-Jf'"*'’^m,n)] 

(l.j)^J’''^  (6) 

Error  Covariance  and  Gain: 

R^"'’''^m,n;k,l)  • I CppR*""  ’ 

(7) 

R^'’'’"^(m,n;m,n)  • 7 C|^^R^'''’’’^(m,n:m-k,n-l)+ 

, (8) 

wheru  ^ is  the  support  of  the  state  vector  s.. 

R^'"’'’ii.J;k,1)  • R^"'*'’^1,J;k,l)  - K^"’*"^m-1, 

n-j).R^">.n)(„,n;k,i),(i,j)  (k,l) 

(9) 

K^'"’"’(1.J)  • R^"'*''hm,n;i,J)/[R^'"*"J(m,n;m,n) 

<1. 


Further  reduction  in  computation  can  be  obtained 
by  computing  (7)  and  (9)  in  a fixed  size  region. 


smaller  than^^’”^.  Such  a region  will  be  refer- 
red to  asT^""’"^  (see  Figure  1). 
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Identification  Algorithm 

The  implementation  of  the  reduced  update  Kalman 
filter  algorithm  (5-10)  requires  the  knowledge  of  the 
Ckj^'s  in  (5).  In  the  following,  an  algorithm  is 
described  for  identification  of  these  unknown  para- 
meters. Let  us  order  all  a column 

vector  c.  Equaction  (1)  can  thCT  be  written  as: 

s(m,n)  • c\,(m,n)  + w(m,n),  (11) 

where  s,  ,U  the  portion  of  the  state  vector  ^ in  . 
the  moMis  active  memory,  the  pseudo-state  vector  . 

Images  are,  in  general,  non- homogeneous  and 
hence  the  elements  of  £ are  spatially  variant. 
Ideally,  then,  a c vector  should  be  found  for  each 
pixel  (m,n).  This,  however,  would  Involve  a large 
amount  of  computation  in  both  the  identification 
and  the  reduced  update  Kalman  filter  algorithms. 

As  a compromise  of  accuracy  versus  amount  of  compu- 
tation, we  assume  regional  homogeneity  and  find  a 
estimate  £ for  each  such  region.  The  least-squares 
estimate  of  c,  for  a KxK  block  of  the  NxN  image,  is 
then  obtaine?  by  minimization  of  the  following  cost 
function:  K,K  ..t  , 

0 . I Cs(1,j)  - c s^(i,j)]‘=  (12) 

K 1 , j 

which  yields  {cf.[3,4]): 

C - [Ps,(1,j)  s^^(1,j)r’-  [^'^sfi,j)s(I.J)] 
iij  ifj^ 


(13) 


Effectively,  we  are  Implementing  a fixed  memory  .vj'Jt* 
filter  over  the  KxK  block  of  an  image. 


The  variance  of  the  plant  noise  w(K)  associated 
with  each  c can  be  est1mated_using: 

°w(K)-^J 


(14) 
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Xhe  parameter  estimate  c and  plant  noise  variance 
o 2 jre  then  used  by  the  reduced  updated  Kalman 

filter  algorlthm.as  explained  In  the  next  section. 
Combined  Identification  and  Estimation 


the  consistency  of  the  Identification  procedure. 
The  test  data  field  was  subsequently  Immersed  In 
noise  and  filtered  with  the  2-0  reduced  update 
filter.  The  computed  and  measured  mean  square 
error  were  within  lOX  for  the  SNR  « 3db  test  case. 


As  previously  mentioned,  the  image  data  field. 

In  general.  Is  not  a stationary  field  and  If  repres- 
sented  by  a constant  set  of  parameters,  a high 
plant  noise  can  be  expected.  On  the  other  hand, 
smaller  blocks  of  the  Image  data  field  are  more 
likely  to  possess  the  statlonarlty  property,  and 
thus  can  be  modelled  with  constant  parameters  more 
accurately. 

Hence  for  simulation  purposes,  the  128  x 128 
image  data  field  was  divided  into  16  32  x 32  blocks. 
Then,  the  1024  values  of  observations  from  each 
block  were  used  to  estimate  the  parameters  of  the 
state  model  parameter  vector  through  (13),  as  well 
as  to  estimate  the  associated  plant  noise  In  each 
block  (14). 

These  parameter  values  and  associated  plant 
noise  were  then  passed  on  to  the  reduced  update 
Kalman  filter  (5»10).  From  previous  experiments, 
the  required'/'^  and(J(^^  regions  were  chosen  as 
shown  in  Fig.  T.  8oundfry  conditions  for  eqns.  (7), 
(8)  and  (9)  were  assumed  to  be  white  Gaussian  while 
those  for  eqn.  (5)  and  (6)  were  assumed  to  be  zero. 
Though  the  parameter  and  gain  values  were  changed 
across  the  boundary  of  each  block  In  the  filter,  no 
detectable  edge  effects  were  noted  In  the  estimated 
Image. 

The  estimated  Image  was  then  compared  with  the 
noise-free  data.  Experimental  results  are  given  In 
the  next  section. 


1 

2 

3 

4 

5 

6 

7 

8 

9 

K) 

11 

12 

• 

• 

iie,ii 


Figure  2:  State  and  Covariance  Update  Regions 
Used  for  Experiments 

Experimental  Results 

To  examine  the  behavior  of  the  Identification 
algorithm.  It  was  first  applied  to  an  Image  gener- 
ated with  a prior  known  parameters,  to  be  explained 
In  part  A of  this  section.  In  part  8,  the  applica- 
tion of  the  Identification  algorithm  to  a typical 
non-homogeneous  Image  Is  discussed.  This  Image, 
observed  through  equation  (4),  was  then  filtered  by 
using  the  reduced  update  Kalman  filter  algorithm 
using  the  Identified  parameters. 

A.  TMtInq  the  Identification  Algorltlwi 

A homogeneous  data  field  was  generated  with 
known  parameters  and  plant  noise  variance,  so  as  to 
satisfy  ni)  for  all  (m,n).  The  results  of  apply- 
ing the  Identification  algorithm  are  given  In  Table 
I.  As  evident  from  these  results,  all  estimates 
Improve  with  an  Increase  In  block  size,  validating 


8.  Processing  of  Noisy  Pictures 

The  noise-free  Image,  as  shown  In  Figure  3,  Is 
a 128  X 128  Image  data  field.  Additive  measure- 
ment noise  w(n)  was  simulated  using  a Gaussian 
white  noise  generating  subroutine.  For  simulation 
purposes,  a 3 db  signal  to  noise  ratio  was  used, 
with  signal  variance  equal  to  2884.  This  nolsy- 
Image  was  processed  as  described  above.  The  proc- 
ess noise  variance  was  experlmentlally  optimized 
to  produce  the  best  MSE  Improvement. 


Figure  4 shows  the  3db  noisy  Image.  The  esti- 
mated Image,  using  a stationary  second  order  ff 
NSHP  models  Is  shown  In  Figure  5.  The  MSE  was 
164  equivalent  to  a 9.44  db  Improvement.  Figure 
6 shows  the  estimated  Image,  using  block  processing 
as  described  In  section  4.  The  MSE  was  134 
equivalent  to  a 10.32  db  Improvement.  This  re- 
presents about  a 20X  reduction  In  MSE.  Subject- 
ively the  noise  level  seems  greatly  reduced,  how- 
ever It  appears  that  some  distortion  has  been  In- 
troduced by  the  filtering  In  some  blocks.  To  In- 
vestigate this  matter,  the  noise  free  Image  and 
the  white  noise  field  were  seperately  filtered 
with  the  results  shown  In  Fig.  7a  and  7b,  respec- 
tively. Figures  7 show  that  the  distortion  Is 
composed  of  colored  noise  that  has  been  shaped  by 
the  recursive  estimator  to  have  a 'local'  spectrum 
similar  to  that  of  the  noise  free  Image.  Such  a 
noise  Is  known  to  be  of  Increased  objectlonality 
over  an  equivalent  amount  of  Isotropic  noise.’ 

This  effect  Is  a direct  result  of  the  filter's 
being  'tuned'  to  the  'local'  spectrum  of  the  signal. 
Thus,  although  subjectively  somewhat  heightened  by 
the  block-wise  constant  model,  this  effect  Is  felt 
to  be  fundamental  to  spatially  varying  filters.  A 
goal  of  future  work  will  be  to  ameliorate  this  ef- 
fect. Possible  approaches  Include  reducing  the 
block  size.  We  also  note  this  effect  Is  most  pro- 
nounced In  flat  portions  of  the  Image. 


These  preliminary  results  are.  In  a sense, 
upper  bounds  on  what  can  be  achieved  In  practice. 
This  Is  because  the  model  parameters  were  calcu- 
lated from  the  noise  free  Image.  For  the  constant 
parameter  model,  the  parameters  could  just  as  well 
be  estimated  from  the  nolsey  Image  (given  the  power 
of  the  white  observation  noise).  However  for  the 
block-wise  constant  parameter  model  with  small 
block  sizes,  use  of  the  nol^  data  would  be  ex- 
pected to  degrade  the  parameter  estimates,  espec- 
ially so  a low  SNR's.  A possible  solution  to  this 
problem  Is  the  use  of  prototype  Images.’  The 
motivation  for  doing  this  Initial  work  Involving 
noise  free  parameter  estimation  was  to  try.  In  so 
far  as  possible,  to  separate  the  effects  of  spati- 
ally varying  filtering  from  the  effects  of  statist- 
ical errors  In  the  parameter  estimates.  Future 
work  will  account  for  the  bias  Introduced  Into  the 
parameter  estimates  by  the  observation  noise  using 
the  method  of  [7]. 
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Conclusions 


The  above  results  show  that  spatially  variant 
recursive  estimation  can  significantly  Improve 
Image  quality  compared  to  spatially  Invariant  or 
constant  parameter  processing.  Conversely, 
spatially  variant  processing  tends  to  correlate 
Input  noise  In  a manner  similar  to  the  'local*  cor. 
relation  of  the  signal,  thus  Increasing  subjective 
objectlonallty  In  flat  portions  of  the  Image. 
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Fig.  5 Estimate  Using  Constant  Model 
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F1g.  7a  Filtered  Noise  Free  Image 


Fig.  7b.  Filtered  White  Noise 
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